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Levy flights in steeper than harmonic potentials have been shown to exhibit finite variance and a 
critical time at which a bifurcation from an initial mono-modal to a terminal bimodal distribution 
occurs (Chechkin et al., Phys. Rev. E 67, 0f0f02(R) (2003)). In this paper, we present a detailed 
study of Levy flights in potentials of the type U(x) oc \x\ c with c > 2. Apart from the bifurcation 
into bimodality, we find the interesting result that for c > 4 a trimodal transient exists due to 
the temporal overlap between the decay of the central peak around the initial ^-condition and the 
building up of the two emerging side-peaks, which are characteristic for the stationary state. Thus, 
for certain system parameters there exists a transient tri-modal distribution of the Levy flight. These 
properties of LFs in external potentials of the power-law type can be represented by certain phase 
diagrams. We also present details about the proof of multi-modality and the numerical procedures 
to establish the probability distribution of the process. 
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I. INTRODUCTION 



Levy flights (LFs) are stochastic, Markov processes, 
which differ from regular Brownian motion by the oc- 
currence of extremely long jumps, whose length is dis- 
tributed according to a Levy stable law with the long 
tail ~ W -1- ", such that its second moment diverges 
ELSIE IS IS- This property strongly contrasts the clas- 
sical Gaussian description of diffusion processes which 
possess finite moments of any given order 0, Q • Given 
their Markov nature, the divergence of the variance of an 
LF might disqualify them as physically meaningful model 
processes for diffusing particles with a finite mass. Yet, 
LFs have important applications to processes, in which 
no finite velocity is required, such as in the energy dif- 
fusion in single molecule spectroscopy . An impressive 
experimental evidence of Levy processes was reported by 
the group of Walther in the study of the position of a 
single ion in a one-dimensional optical lattice, in which 
diverging fluctuations could be observed in the kinetic 
energy [9|- From a phenomenological point of view, LFs 
have been used to describe the dynamics observed in plas- 
mas ^3 or m molecular collisions They have also 
been successfully applied to describe the statistics en- 
countered in the spatial gazing patterns of bacteria |12| , 
albatross birds |l^. or even spidermonkeys Q. Prob- 
ably the earliest application of LFs, however, may be 
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in the modelling of financial markets |15| . LFs were 
shown to give rise to surprisingly rich band structures 
in periodic potentials |lf|. Reverse engineering methods 
have been developed to construct a Langevin system with 
Levy noise to produce a pre-defined steady state 01 • We 
also note that the stationary state of LFs in a confining 
external potential is expected to approach to the one of 
Levy walks for long times (Tsl ITflj (the spatiotemporally 
coupled version of LFs which have converging moments 
of any order) , and therefore the more straightforward de- 
termination of the solution for LFs can be used to gain 
insight into Levy walks. 

The theory of homogeneous LFs is pretty well under- 
stood. Thus, LFs in the continuum limit can, inter alia, 
be described by continuous time random walks with long- 
tailed, asymptotic power-law jump length distributions 
|18| , or alternatively by a Lang evin equation with 5- 
correlated, Levy noise [2(1 lU • Both descriptions can 
be map ped onto the (space-) fractional diffusion equation 



df d a 
1 =D—f(x,t), 



0t 



d\x\ 



(1) 



where the fractional Riesz operator d a /d\x\ a is most eas- 
ily defined in terms of its Fourier transform |27j 



ikx 



f(x,t)dx = -\k\ a f(k,t); 



(2) 



here, < a < 2. In the following we restrict our analysis 
to the range 1 < a < 2. In the limit a = 2, naturally, 
equation Q reduces to the classical diffusion equation, 



2 



which described Gaussian transport. From the Fourier 
transform of equation QJ, df/dt = — \k\ a f(k, t), we con- 
clude that the characteristic function is 



f{k,t) = exp{-Dt\k\ a ) 



(3) 



i.e., exactly the characteristic function of a symmetric 
Levy stable law of index a 0,0]. In position space, ex- 
pression © can be represented exactly in terms of Fox 
H- functions |2£|, but for our purposes in what fol- 
lows, it is enough to remember that for large \x\, we 
find f{x,t) ~ Dt/\x\ 1+a , such that the variance diverges: 
(x 2 (t)) — oo. Also the general formulation of LFs in ex- 
ternal potentials and in phase space seems well founded, 
in terms of fractional Klein-Kramers equations; see, e.g., 
references [H |H HI EJ 

However, much less is known about the actual be- 
haviour of LFs in external potentials U(x), i.e., about 
the properties of the probability density function (PDF) 
f(x,t) in the presence of such a U(x). One of the few 
solved examples is the one of LFs in a harmonic poten- 
tial [33, whose solution always follows the same stable 
law of index a, like in the regular Ornstein-Uhlenbeck 
process the PDF always stays Gaussian. In general, 
LFs in the presence of an external potential field in 
the overdamped case are governed by the (space-) frac- 
tional Fokker-Planck (or Einstein-Smoluchowski) equa- 
tion EJIHHIi! 



df_ _ ( d U'jx) 
dt \ dx mrj 



D 



(4) 



which has the characteristic property that the drift term 
enters with the usual first order derivative and thus pre- 
serves its additive quality such that for a constant force 
F(x) — Vmr) the solution is given by the drift-free PDF 
taken at the similarity variable x — Vt, and for a general 
external potential the stationary solution differs from the 
Boltzmann distribution. This latter property is in con- 
trast to an alternative LF model suggested in reference 

An a priori unexpected behaviour of LFs was obtained 
recently as solution of the fractional Fokker-Planck equa- 
tion namely the occurrence of bimodal solutions 
(PDFs with two maxima) and the finiteness of the sec- 
ond moment in the presence of superharmonic potentials 
of the form SHil 



U(x) 



„2m+2 



2m + 2 



1,2, 



(5) 



Note that the amplitude a > has dimension [a] = 
g ■ cm _2m /sec 2 . Thus, for a potential of the form U (x) = 
ax 2 + bx A with a, b > 0, a turnover can be tuned from 
the properties of the solution of the harmonic problem 
(mono-modal, diverging variance) to a finite variance and 
bimodal solution by varying the ratio b : a |35|. If such 
bimodality occurs, it results from a bifurcation at a crit- 
ical time t c . A typical result is shown in figure ^ 



for the quartic case m = 1 and Levy index a = 1.2: 
from an initial 5-peak, eventually a bimodal distribution 
emerges. Note that we use dimensionless quantities in 
the plots, as introduced below. The location of the global 
maximum/maxima is displayed in figure |2 where the bi- 
furcation is a distinct mark. In the same figure, the value 
of the PDF / at the newly emerging humps is compared 
to the value at the origin. 

In what follows, we briefly review the description of 
LFs in external fields in terms of the Langevin equation 
with white Levy noise and the space-fractional Fokker- 
Planck equation, as well as the properties of LFs in su- 
perharmonic potentials; and we present two proofs, for 
both the finite variance in these cases and the existence of 
multi-modality. We then proceed to show that for poten- 
tials steeper than JSJ with m = 1, even a tri- modal PDF 
can be obtained. In two phase diagrams, we can classify 
the existence of the different n-modal states in depen- 
dence of the potential exponent c and the Levy index a; 
and the critical bifurcations between different n-modal 
domains as a function of time t and potential exponent 
c. In the appendix, we discuss the numerical methods 
from which the PDFs are obtained. 



II. STARTING EQUATIONS 

In this section, we formulate the dynamical description 
of LFs on the stochastic differential (Langevin equation) 
and the deterministic (fractional Fokker-Planck equa- 
tion) levels. For the latter, we also discuss the corre- 
sponding form in Fourier space. 



A. Starting equations in real space 

1. The Langevin equation with Levy noise 

On the level of the stochastic description, our starting 
point is the overdamped Langevin equation [23, 0, |3^| 



dx F{x) 



Y a (t), 



dt M-f 

where F — —dU/dx is an external force with potential 

a\x'" 



(6) 



U{x) 



with a > and c > 2, M is the particle mass, 7 the fric- 
tion coefficient, and Y a (t) represents a stationary white 
Levy noise with Levy index a (1 < a < 2). 

We employ the white Levy noise Y a (t) such that the 
process 



ft+At 

L(At) - / Y a (r)dT, 



(8) 



which is the time integral over an increment At is an te- 
stable process with stationary independent increments. 



3 



Restricting ourselves to symmetric Levy laws, this im- 
plies the characteristic function 



with 



p L (k,At) =cxp{~D\k\ a At). 



(9) 



The constant D in this description has the meaning of the 
intensity of the Langevin source, and [D] = cm™/ sec. 

In figure |3 we show realisations of white Levy noises 
for various values of a. The sharply pronounced 'out- 
liers', due to the long-tailed nature of the Levy stable 
distribution, are distinct, in comparison to the Gaussian 
case a = 2. 



2. Fractional Fokker-Planck equation 

The Langevin equation © is still of the Markov- 
type, and it is therefore fairly straightforward to show 
that the corresponding fluctuation-averaged (determinis- 
tic) description is given in terms of the space-fractional 
Fokker-Planck equation © [U EU, El . In what follows, 
we derive solution of equation Q for the J-initial condi- 
tion 



f(x,0) = 6(x). 



(10) 



The space-fractional derivative d a /d\x\ a occurring in 
the fractional Fokker-Planck equation (jJJ is called the 
Riesz fractional derivative, defined through 



d a f 
d\x\ a 



2 cos(7ra/2) ' / 

iHf, a=l 



where we use the following abbreviations: 



and 



(D°/)(a:) 



r(2 - a) dx 2 I M (a; - ' 



r(2 - a) dx 2 J x {£ - x)"- 1 



(11) 



(12) 



(13) 



for, respectively, the left and right Riemann-Liouville 
derivatives (1 < a < 2); and [3(| 



(H/)(x) 



1 



x-a 



(14) 



is the Gilbert transform operator. The definitions for 
d a /d\x\ a demonstrate the strongly non-local property of 
the space- fractional Fokker-Planck equation, i.e., strong 
correlations in x. 



3. Rescaling of the dynamical equations 



Passing to dimensionless variables 

x' = x/xq, t = t/to, 



(15) 



(MD~{\ 



l/(c~2+a) 



«o=^j , t = ^, (16) 

the starting equations take the form (we omit primes 
below) 



dx dU T , , . 
-dt--^ + Y ^ 

instead of the Langevin equation 10, and 

df(x,t) d <iu d a f 



(17) 



(18) 



dt dx dx d\x\ a 
instead of the fractional Fokker-Planck equation and 



U(x) = L-L, 
c 



B. Starting equations in Fourier space 



(19) 



If f(k,t) denotes the characteristic function (CF), i.e., 
Fourier transform of f(x,t), we write 



f(x,t) + f(k,t), 



(20) 



where we use the sign -j- to denote a Fourier transform 
pair. Since 27] 



and 



we obtain 



{T>%f){x,t) + {Tik) a f{k,t), 



(H/)(z,i)-isign(AO/(M), 



d a f 
d\x\ a 



\k\ a f(k,t), 



(21) 



(22) 



(23) 



for all a's. The equivalent of the fractional Fokker-Planck 
equation H18|) for the CF then follows immediately, 



dt 



\k\ a f = U fe J, 



with the initial condition 

f(k,t = 0) = 1, 
and the normalisation 

/(fc = 0,t) = l. 



(24) 

(25) 
(26) 



The external potential U (x) turns into the linear differ- 
ential operator in k, 



U fc / = 



d fdU 



dx V dx 



ik I e ikx sign(x)\x\ c - L f(x,t)dx. (27) 



f)dx 
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Next, making use the following 'inverse' expressions 

<±ix) a f(x) + (D%f)(k), (28) 

and 

-i(sign(x)f(x) + (Hf)(k). (29) 

we obtain the explicit expression for the external poten- 
tial operator, 



U fc / 



2cos(ttc/2) 



(^-■D- )/, c ^3, 5, 7,. 



{-l) m k4^Hf 



c = 3, 5, 7, . 



(30) 

Note that for the even potential exponents c = 2m + 2 , 
rn = 0, 1, 2, . . ., we find the simplified expression 



2m+l 



U fc = (-l) m+i fc 



dk 2m+1 ' 



(31) 



We see that the force term can be written in terms of 
fractional derivatives in Fourier space, and therefore it is 
not straightforward to calculate even the stationary solu- 
tion of the fractional Fokker-Planck equation l|18[) in the 
general case c ^ N. In particular, in this latter case, the 
non-local equation (|18|l in x-space translates into a non- 
local equation in fc-space, where the non-locality shifts 
from the diffusion to the drift term. 



III. ANALYTICAL RESULTS 

In the preceding section, we discussed some elementary 
properties of the space-fractional Fokker-Planck equation 
for LFs, in particular, we pointed out the spatially non- 
local character of equation l|18l) . and its Fourier space 
counterpart l|24|) . In this section, we determine the ana- 
lytical solution of the fractional Fokker-Planck equation. 
We start with the exactly solvable stationary quartic 
Cauchy oscillator, to demonstrate directly the occurring 
bimodality, and then move on to the general case. The 
major results will be the determination of rt-modality, 
finite variance, and the parametric dependence of the as- 
sociated bifurcations. 



in Fourier space. Its solution is 

whose inverse Fourier transform results in the simple an- 
alytical form 



1 



7r(l - X 2 + X 4 ) ' 



(35) 



As shown in figure^] this solution has two global maxima 
at a; max = ±l/\/2 apart from a local minimum at the 
origin (the position of the initial condition, that is) and 
its variance 

(x 2 ) = 1 (36) 

is finite, due to the long-tail asymptotics f s t(x) ~ x~*. 
These two distinct properties for LFs will turn out to be 
a central theme in the remainder of this work. 



B. Formal solution of equation ()18[l 

Returning to the general case, we rewrite equation l|24|) 
in the equivalent integral form, 



f(k,t)=p a (k,t) + dTPa(k,t-T)Vkf(k,T) (37) 

Jo 



where 



p a (k,t) = exp(-|fc| a i) 



(38) 



is the CF of a free (homogeneous) LF. This relation fol- 
lows from equation (|24|) via formally treating it as a non- 
homogeneous linear differential equation of first order, 
where Ufe plays the role of the non- homogeneity. Then, 
()24J) is obtained from variation of constants. (Differenti- 
ate equation l|37|) to return to (|24|l .') 

Equation (|37|l can be solved formally by iterations: Let 



/ (0) (M)=£a(M), 



(39) 



then 



A. The stationary quartic Cauchy oscillator 

Let us first regard the case of the stationary quartic 
potential with c = 4 for the Cauchy-LF with a = 1, i.e., 
the solution of the equation 



d x 3 f st {x) + -S-rUix) = 0, 
a\x\ 



dx 



d 3 f st (k) 
dk 3 



sign(fe)|fc|/ st (fe) 



(32) 



(33) 



/ (1) (fc,<) =p a (k,t)+ / dTp a (k,t-T)U k f (0) (k,T), 

Jo 

(40) 

/ (2) (M) =P Q (fc,i)+ / d T p a {k,t~T)\J kPa {k,T) + 



dr / dT'p a (k,t - T)\J kPa (k,T - T')U k p a (k,T')(41) 
o Jo 



etc. Invoking the definition of the convolution, 

A*B= [ drA(t — t)B(t) = [ dTA{r)B{t- r), (42) 
Jo Jo 
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and using 

A*B*C = (A*B)*C = A*(B* C), (43) 
we arrive at the formal solution 



f(k,t) = ^Tp a (*u k p a ) n , 



(44) 



This procedure is analogous to perturbation theory, Ufc/ 
playing the role of the interaction term, see, for instance, 
reference j33, chapter 16. 

Applying a Laplace transformation, 



/(M= / dt exp(-st)/(M), (45) 
Jo 

to equation (|44(1 . we obtain 

}(k, s) = f a (k, s) + f a (k, s)V k f(k, s), (46) 



where 



k a 



(47) 



is the Fourier-Laplace transform of the homogeneous 
Levy stable PDF. Thus, we obtain the equivalent of so- 
lution (|4*4*jl in (fc, s)-space: 



oo n 

f(k,s) = ^2 [p*( k ,s)Uk P a {k,s) 



n=0 



(48) 



This iterative construction scheme for the solution of the 
fractional Fokker-Planck equation will turn out to be use- 
ful below. 



C. Existence of a bifurcation time. 

For the case of the unimodal initial condition f(x, 0) = 
6(x) we now proof the existence of a finite bifurcation 
time ti2 for the turnover from unimodal to bimodal PDF. 
At this time, the curvature at the origin will vanish, i.e., 
be an inflection point: 



5V 
dx 2 







X—0.t — ti2 



Introducing 



J(t) = / dkk 2 f(k,t), 



(49) 



(50) 



equation l|49|) is equivalent to (note that the CF is an 
even function) 



o. 



(51) 



The bifurcation can now be obtained from the iterative 
solution l|48l) : we consider the specific case c = 4. From 
the first order approximation 



and we have 



s + k c 



U k = k 



s + k° 



dk 3 ' 



(52) 



(53) 



Combining these two expressions, we produce 

~ \ k a ~ 2 
(fc, s) = ; h a(a — 1)(2 — a)- — -^r 



+6a 2 (a - 1)- 



k 2 



6a d 



k 



3a-2 



(s + k a ) 4 (s + k a ) 5 

or, after inverse Laplace transformation, 



. (54) 



A(M)=e 



-k"t 



1 - —t 4 k 3a - 2 + a 2 (a - l)t 3 k 2a - 2 
4 



+a(a- 1){2- a) — k a ~ 2 



(55) 



The first approximation to the bifurcation time ti2 is 
then determined via equation (|50|l . i.e., we calculate 



to obtain 



' 12 



dk k 2 f l {k,t { l 2 ) )=Q, 



4T(3/a) \ Q/(2+Q) 



(56) 



(57) 



v 3(3-a)r(l/a). 

In figure El we show the dependence of this first approx- 
imation as a function of the Levy index a (dashed 
line), in comparison to the values determined from the 
numerical solution of the fractional Fokker-Planck equa- 
tion (|18|1 shown as the dotted line. The second order iter- 
ation for the PDF, f2(k, t), can be obtained with maple6, 
from which in turn the second approximation for the bi- 
furcation time is found in analogy to above procedure. 
The result is displayed as the full line in figure |SJ The 
two approximative results show in fact surprisingly good 
agreement with the numerical result of the full PDF. 



D. Proof of non-uni-modality of stationary solution 

for c> 2 

In this subsection we demonstrate that the stationary 
solution of the kinetic equation l|18|l has a non-unimodal 
shape. For this purpose, we use an alternative expression 
for the fractional Riesz derivative (compare, e.g., refer- 
ence 



d a f(x) 
d\x\ a 



r(l + a) 



sin(a7r/2) 



„ f(x + Q-2f(x) + f(x-Z) 
dt, (58) 
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valid for < a < 2. This definition can be obtained by 
'regularisation' of our definitions to for a / 1; 
however, expression H58JI is also valid at a = 1. In the 
stationary state (df /dt — 0), we get from equation (|18|l : 



^ (8g „<*r'/..) + ^f = o. 

Thus, it follows that at c > 2 (strict inequality) 



d\x\ 



0, 



(59) 



(60) 



x=0 



or, from definition (|58fl and taking into account that 
fst{x) is an even function, 



/ st (0-/st(o) 



l + Q 



0. 



(61) 



From this latter relation, we can immediately obtain 
proof the non-uni- modality of / s t, which we produce in 
two steps: 

(1) If we assume that the stationary PDF f s t(x) is 
unimodal, then due to the symmetry x — > —x, it neces- 
sarily has one global maximum at x = 0. In this case 
the integrand in equation l|61[) must be negative, and 
therefore contradicts equation If) . Therefore, f s t(x) is 
non-unimodal. 

(2) We can in addition exclude /(0) = 0, as in this case 
the integrand will be positive, which is again in contra- 
diction with equation l|61l) . 

Since f(x) — > at x — > oo, basing on statements (1) 
and (2), one may conclude that the simplest situation is 
such that £o > exists with the property 



and, 



/(0) 



>o, 



(62) 



(63) 



i.e., the condition for two- hump stationary PDF for all 
c > 2. At intermittent times, however, we will show that 
also a tri- modal state may exist. 



E. Power-law asymptotics of stationary solutions 
for c > 2, and finite variance for c > 2 

We now derive the power-law asymptotics of the sta- 
tionary PDF f s t(x) for external potentials of the form 
(JSJ) with general c > 2. To this end, we note that at 
x — > +oo, it is reasonable to assume 



D a _f st « D°f m 



(64) 



since the region of integration for the right-side Riemann- 
Liouville derivative (D"/ St )(2:), (x,oo), is much smaller 



than the region of integration for the left-side derivative 
(D"/ St )(a;), (— oo, x), in which also the major portion of 
fst(x) is located. Thus, at large x we get for the station- 
ary state, 



d 

dx 



dU 
dx 



1 



M ) 2cos( 7ra /2)<W- co (z-O"- 1 ' 

(65) 

Equivalently, this corresponds to the approximative 
equality 



1 



2 cos(7ra/2) dx 



z-£) a - ] 



(66) 



We are seeking asymptotic behaviours of f s t(x) in the 
form f{x) w C\/x^ (x — > +oo, /i > 0). After integration 
of relation l|66|) . we find 



2CiCOs(7ra/2)r(2-a) 
—fJ. + c 



x-£) a - 



(67) 



The integral on the right hand side can be approximated 
through 



1 



/ st (£R = 



l 



f 



(68) 



Thus, we identify the powers of x and the prefactor, with 
the results 



(j, = a + c — 1 



and 



sin(7ra/2)r(a) 



(69) 



(70) 



By symmetry of the PDF we therefore recover the general 
asymptotic form 



fix) 



sin(7ra/2)r(a) 
7r|aj| M 



+oo 



(71) 



for all c > 2. This result is remarkable, for various rea- 
sons: 

(i) despite the approximations involved, the asymp- 
totic form l|71(l for arbitrary c > 2 matches exactly previ- 
ously obtained forms, such as the exact analytical result 
for the harmonic LF (linear Levy oscillator), c = 2 re- 
ported in reference [32j; the result for the quartic Levy 
oscillator with c = 4 discussed in references [3 l35| : 
and the case of even power-law exponents c = 2m + 2 
(to e No) given in reference [34| . 

(ii) The prefactor C\ is independent of the potential 
exponent c; in this sense, C\ is 'universal'. 

(iii) For each value a of the Levy index the 'critical' 
value 

c cr = 4 - a (72) 

exists such that at c < c cr the variance < x 2 > is infinite, 
whereas at c > c cr the variance is finite. 

(iv) We have found a fairly simple trick to construct 
stationary solutions at large x in the form of inverse 
power series. 
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IV. NUMERICAL RESULTS 

In this section, we show results for the PDF of the 
fractional Fokker-Planck equation (|18fl obtained via two 
different numerical techniques, one being the Grunwald- 
Lctnikov method, which is based on an iterative solution 
of the deterministic dynamical equation (|18|) by replacing 
the Riesz fractional derivative with Griinwald-Letnikov 
operators; the other being the Langevin method, in which 
the Langevin equation with Levy noise, equation J6j is 
integrated numerically. Both methods lead to analogous 
results, and they also produce results for the PDF which 
are in perfect agreement with above analytical results. 
We present the numerical results in two subsections, de- 
voted to the two numerical methods. These methods 
themselves are discussed in the appendix. 



A. Results from the Griinwald-Letnikov method 

1. Tri-modal transient state at c > 4 

Before, we have proved the existence of a bimodal sta- 
tionary state for the quartic (c = 4) Levy oscillator. This 
bimodality emerges as a bifurcation at a critical time 
t CT , at which the curvature at the origin vanishes. This 
scenario is replaced for c > 4, as displayed in figure 
Thus, there obviously exist two timescales, the critical 
time for the emergence of the two off-centre maxima, 
which are characteristic for the stationary state; and a 
second one, which is designated to the relaxation of the 
initial central hump, i.e., the decaying initial distribu- 
tion f(x,0) = S(x). The formation of the two off-centre 
humps while the central one is still present, is detailed 
in figure [7| This existence of a transient tri-modal state 
was found to be typical for all c > 4. 

In figures and |§1 we show additional details of the tri- 
modal state. Thus, figure|S|depicts a bifurcation diagram 
for the abovely shown process; the initial mono-modal 
PDF bifurcates to a tri-modal one, before it finally turns 
over to bimodality. In figure El these two turnover times 
are displayed as function of the Levy index a. Clearly, 
there is always a gap between these two time scales, leav- 
ing the intermittent time for the tri-modal state, and for 
a — > 2, this region shrinks, both curves converge to infin- 
ity, as in the regular Gaussian case no such bifurcation 
exists. 



2. Phase diagrams for n-modal states 

Above findings can be put in context with the purely 
bimodal case discussed earlier. A convenient way of dis- 
playing the n-modal character of the PDF in the presence 
of a superharmonic external potential of the type (|19|l is 
the phase diagram shown in figure 1101 There, we sum- 
marise the findings that for 2 < c < 4 the bifurcation 



occurs between the initial mono-modal and the station- 
ary bimodal PDF at a finite critical time, whereas for 
c > 4, a transient tri-modal state exists. Moreover, we 
also include the shaded region, in which c is too small to 
ensure a finite variance. Complementarily, in figure 1111 
the temporal domains of the n-modal states are graphed, 
and the solid lines separating these domains correspond 
to the critical timescales. Again, the transient nature of 
the tri-modal state is distinct. 



B. Langevin method 

As explained in the appendix, this method directly in- 
tegrates the Langevin equation for white Levy noise, the 
latter being portrayed in figure 01 Typical results for the 
sample paths under the influence of an external potential 
(|19H with increasing superharmonicity are shown in figure 
I12l in comparison to the Brownian case (i.e., white Gaus- 
sian noise). For increasing the external exponent c, the 
long excursions, which are typical for homogeneous LFs 
are increasingly suppressed. In the harmonic case c = 2 
still present (in this case, the variance is diverging), they 
are clearly confined for c > 2. Note the comparable ordi- 
nate windows in comparison to the significantly different 
scale in the homogeneous case of figure |3] For all dis- 
played cases, however, the qualitative behaviour of the 
noise under the external potential is different from the 
Brownian noise even in this case of strong confinement. 
In the same figure, we also show the curvature of the 
external potential. Additional investigations have shown 
that the maximum curvature is always very close to the 
positions of the two maxima, leading us to conjecture 
that they are in fact identical. 

The latter observation is further investigated in figure 
1131 On a linear scale, the potential well and its curvature 
are compared to the stationary PDF, clearly demonstrat- 
ing the proximity of maximum curvature and the two 
maxima. Figure El a ls° corroborates on the basis of the 
Langevin method the asymptotic inverse power-law be- 
haviour derived in equation (|71() . 

Finally, in figure El we display the time evolution of 
the PDF in the three different modality-regimes accord- 
ing to figure im The comparatively noisy result is due to 
a small number of trajectories used for the statistical av- 
erage, due to the rather computation intensive program. 



V. CONCLUSIONS 

By combining analytical and numerical results, we 
discussed LFs in a superharmonic external potential of 
power c. Depending on the magnitude of this exponent 
c, different regimes could be demonstrated. Thus, for 
c = 2, the character of the Levy noise imprinted on the 
process, is not changed by the external potential: the 
resulting PDF has Levy index a, the same as the noise 
itself, and will thus give rise to a diverging variance at 
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all times. Conversely, for c > 2, the variance becomes 
finite if only c > c cr = 4 — a. This is due to the fact 
that the PDF leaves the class of Levy stable PDFs and 
acquires an inverse power-law asymptotic behaviour with 
power /i = a + c — 1 . Obviously, moments of higher order 
still diverge. Apart from the finite variance, the PDF is 
distinguished by the observation that it bifurcates from 
the initial mono-modal to a stationary bimodal state. If 
c > 4, there exists a transient tri-modal state. This rich- 
ness of the PDF both during relaxation and at station- 
arity, depending on a competition between Levy noise 
and steepness of the potential contrasts the universal ap- 
proach to the Boltzmann equilibrium, solely defined by 
the external potential, in classical diffusion. 

It may be speculated what the exact kinetic reason for 
the occurrence of the multiple humps is. Due to the ob- 
servation that the non-transient humps seem to coincide 
with the positions of maximum curvature of the external 
potential, one may conclude that the potential at these 
points changes almost abruptly (especially for larger c) 
from a rather flat to a very steep slope, towards which 
the random walker is driven by the diffusivity. I.e., in 
comparison to the harmonic and subharmonic cases, the 
restoring force close to the origin is comparably weak. 

The different regimes for c > 2 can be classified in 
terms of critical quantities, in particular, the bifurcation 
time(s) t CT and the critical external potential exponent 
c cr . LFs in super harmonic potentials can then be conve- 
niently represented by phase diagrams on the (c, a) and 
(c, t cr ) plains. 

The numerical solution of both the fractional Fokker- 
Planck equation in terms of the Griinwald-Letnikov 
scheme used to find a discretised approximation of the 
fractional Riesz operator shows reliable convergence, as 
corroborated by direct solution of the corresponding 
Langevin equation. 

Our findings have underlined the statement that the 
properties of LFs, in particular under non-trivial bound- 
ary conditions or in an external potential are not fully 
understood. The general difficulty, which hampers a 
similarly straightforward investigation like in the regular 
Gaussian or the subdiffusive cases, is connected with the 
strong spatial correlations of the problems, manifested in 
the integrodifferential nature of the Riesz fractional oper- 
ator. For this reason it is already non-trivial to determine 
the stationary solution of the process. We expect, by the 
fact that diverging fluctuations appear to be relevant in 
physical systems, a range of yet unknown properties of 
LFs remain to be discovered. 
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Appendix A: NUMERICAL SOLUTION 
METHODS 



In this appendix, we briefly review the numerical tech- 
niques, which were used in this work to determine the 
PDF from the fractional Fokker-Planck equation (|f 8fl and 
the Langevin equation (|f 7(1 . 



1. Numerical solution of the fractional 
Fokker-Planck equation (1181) via the 
Griinwald-Letnikov method 

From a mathematical point of view, the fractional 
Fokker-Planck equation l|f 81) is an first-order partial dif- 
ferential equation in time, and of non-local, integro- 
differential kind in the position co-ordinate x. It can be 
solved numerically via an efficient discretisation scheme 
following Griinwald and Letnikov [3^ . Eol l4l| . 

Let us designate the force component on the right hand 
side of equation (|f 8|l as 



and the diffusion part as 



D(x,t) 



d a f 
d\x\ a 



(Al) 



(A2) 



With these definitions, we can rewrite equation (|f 8fl in 
terms of a discretisation scheme as 



fj,n+l fj,n 

At 



where we encounter the term 



F - x c ~ 2 

r J,n — Xj 



(c- l)fj <n +Xj 



/j+l,7i Jj — l,n 

2Ax 



(A3) 



(A4) 



which is the force component of the potential U(x) — 
\x\ c /c. Here, At and Ax are the finite increments in 
time and position, such that i„ = ndt and Xj — jAx, for 
n = 0, 1, . . . , N and j = 0, 1, . . . , J, and fj >n = f(xj,t n ). 
Due to the inversion symmetry of the kinetic equation 
(|f 8fl . it is sufficient to solve it on the right semi-axis. In 
the evaluation of the numerical scheme, we define xj such 
that the PDF in the stationary state is sufficiently small, 
say, 1 0~ 3 , as determined from the asymptotic form f7T)l . 

In order to find a discrete time and position expres- 
sion for the fractional Riesz derivative in equation ( |A2I) . 
we employ the Griinwald-Letnikov scheme |39l rUi l4lj , 
according to which we obtain 
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2(Aa;) Q cos(7ra/2) 



(A5) 
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where 



with 



1, 



l)...(a-g + l)/g! 



g > 
q < 



(A6) 



(A7) 



and 1 < a < 2. Note that in the limiting case a = 2 only 
three coefficients differ from zero, namely, £o = 1, £i = 
—2, and £2 = 1, corresponding to the standard three- 
point difference-scheme for the second order derivative, 
d 2 g(xj)/dx 2 w - 2g 3 + g j _ l )/(Ax) 2 . In figureCHl 

we demonstrate that with decreasing a, an increasing 
number of coefficients contribute significantly to the sum 
in equation l|A5jl . This becomes particularly clear in the 
logarithmic representation in the bottom plot of figure 
ITol We note that the condition 



fi=At/(Ax) a < 0.5 



(A8) 



is needed to ensure the numerical stability of the dis- 
cretisation scheme. In our numerical evaluation, we use 
Ax = 10~ 3 , and therefore the associated time increment 
At ~ 1CP 5 . . . 10~ 6 , depending on the actual value of a. 
The initial condition for equation (| A3|> is /o,o = 1/ Ax. 

In figureEl the time evolution of the PDF is shown to- 
gether with the evolution of the force and diffusion com- 
ponents defined by equations (jAlfl and l|A2(l , respectively. 



Accordingly, at the initial stage of the relaxation pro- 
cess, the diffusion component prevails. The force term 
grows in the course of time, until at the stationary state 
F — > —D. This is particularly visible at the bottom right 
part of figure [TBI which corresponds to the stationary bi- 
modal state shown to the left. 



2. Numerical solution of the Langevin equation ||{>J| 



An alternative way to obtain the PDF is to sample the 
trajectories determined by the Langevin equation ©. To 
this end, equation l|17|) is discretised in time according to 



X n +1 



F(x n )At + (At) 1/a Y a (nAt) 



(A9) 



with t n — nAt for n = 0,1,2,..., and where the 
force F(x n ) is the dimensionless force field at position 
x n . The sequence {Y a (nAt)} is a discrete-time approx- 
imation of a white Levy noise of index a with a unit 
scale parameter. That is, the sequence of independent 
random variables possessing the characteristic function 
p = exp(— \k\ a ). To generate this sequence {Y a (nAt)} : 
we used the method outlined in reference |42j . 
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Figure 1: Time evolution of the LF-PDF in the presence of 
the superharmonic external potential JSJ with m = 1 (quar- 
tic Levy oscillator) and Levy index a = 1.2, obtained from 
the numerical solution of the fractional Fokker-Planck equa- 
tion, using the Griinwald-Letnikov representation of the frac- 
tional Riesz derivative (full line). The initial condition is a 
5-function at the origin. The dashed lines indicate the corre- 
sponding Boltzmann distribution. The transition from one to 
two maxima is clearly seen. This picture of the time evolu- 
tion is typical for 2 < c < 4 where c is defined in equation 
below. The corresponding location of the maximum/maxima 
as a function of time is shown in figure |2] 



Figure 5: Bifurcation time t\2 versus Levy exponent a 
at external potential exponent c = 4.0. Black dots: bi- 
furcation time deduced from the numerical solution of the 
fractional Fokker-Planck equation 1181 using the Griinwald- 
Letnikov representation of the fractional Riesz derivative (see 
appendix). Dashed line: first approximation solid line: 
second approximation tjj . 



Figure 2: Bifurcation diagrams for the case m = 1.0, a = 
1.2, corresponding to the PDF shown in figure Left: the 
thick lines show the location of the maximum, which at the 
bifurcation time ti2 = 0.83 ± 0.01 turns into two maxima. 
Right: the value of the PDF in the maxima location (thick 
line) and the value in the minimum at x = (thin line). 

Figure 3: White Levy noises with the Levy indexes a = 
2, 1.7, 1.3, 1.0. The 'outliers' are increasingly more pro- 
nounced the smaller the Levy index a becomes. Note the 
different scales on the ordinates. 



Figure 6: Time evolution of the PDF governed by the frac- 
tional Fokker-Planck equation 1181 in a superharmonic po- 
tential @ with exponent c = 5.5, and for Levy index a — 
1.2; obtained from numerical solution using the Griinwald- 
Letnikov method explained in the appendix. Initial condition 
is f(x, 0) = 8(x). The dashed lines indicate the corresponding 
Boltzmann distribution. The transitions between 1 — > 3 — > 2 
humps are clearly seen. This picture of time evolution is typ- 
ical for c > 4. On a finer scale, we depict the transient tri- 
modal state in figure |7| 



Figure 4: Stationary PDF 1351 of the Cauchy-LF in a quar- 
tic (c = 4) potential. Two global maxima exist at x max = 
±^/l/2, and there is a local minimum at the origin. 



Figure 7: The transition 1 — » 3 — » 2 from figure [(Jon a finer 
scale (c = 5.5, a = 1.2). 



Figure 8: Bifurcation diagrams for the case c = 5.5 and 
a = 1.2 corresponding to figures E] and Q Left: positions 
i mx of the maxima (global and local, thick lines); the thin 
lines indicate the positions of the minima (at the first bifur- 
cation time, there is a horizontal tangent at the site of the 
two emerging off-centre maxima). The bifurcation times are 
tia = 0.75 ± 0.01 and t 32 = 0.92 ± 0.01. Right: values of the 
PDF at the maxima (thick lines); the thin line indicates the 
value of the PDF in the minima. 

Figure 9: Bifurcation times ii3 versus a (lower curve) and t%2 
(upper curve) for the external potential exponent c = 5.5. 



Figure 10: (c, «)-map showing different regimes of evolution 
of the PDF, and the stationary states. The region with infi- 
nite variance is shaded. The region c < 4 covers the transition 
from 1 to 2 humps during time evolution. For c > 4, a tran- 
sition from 1 to 3, and then from 3 to 2 humps occur. For all 
c's there are two maxima in the stationary state. Compare 
figure 1111 

Figure 11: (c, t)-map showing states with different number 
of maxima, and the transitions between them during time 
evolution. Region 1: PDF has 1 hump; region 2: PDF has 2 
humps; region 3: PDF has 3 humps. At c < 4 there is only a 
transition 1 — » 2, whereas at c > 4 there are two transitions: 
1 — > 3 and after the transient tri-modal regime, 3^2. 



Figure 12: Left column: the potential energy functions 
U = x c /c, (solid lines) and their curvatures (dotted lines) for 
different values of c: c = 2 (linear oscillator), and c = 4, 6, 8 
(strongly non-linear oscillators) . Middle column: typical sam- 
ple paths of the Brownian oscillators, a = 2, with the poten- 
tial energy functions shown on the left. Right column: typical 
sample paths of the Levy oscillators, a — 1. It is seen that 
with c increasing the potential walls become steeper, and the 
Levy flights become shorter; in this sense, they are "confined" . 



Figure 13: Stationary PDF f B t(x) on linear (left) and double- 
logarithmic (right) scale, obtained from the Langevin equa- 
tion for a) c = 4.0, b) c = 5.7, and c) c = 6.5. The thin 
lines on the left show the potential wells and their curvatures. 
The solid lines on the right show the asymptotics as given by 
equation I71H . 



Figure 14: Time evolution of the PDF f(x) obtained from the 
Langevin equation for c = 5.5, a = 1.2; besides the statistical 
averaging over different trajectories, a time averaging over a 
small time interval has been used to create these figures: a) 
inside the region of a single peak, t < tis; b) inside the region 
of three peaks, tis < t < £32; c) inside the region of two peaks, 

t > t 32 - 



13 



Figure 15: Coefficients £ 9 in Griinwald-Letnikov approxima- 
tion for different values of the Levy index a = 1.9, 1.5, and 
1.1. 



Figure 16: Further details of the Griinwald-Letnikov scheme. 
Left: Time evolution of the PDF as obtained by numerical 
solution of equation (IA3II at c = 4 and a — 1.2. Right: Time 
evolution of the diffusion component 1A2I (thick lines), and 
the force term (IA1I> (thin lines). 
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